Direct Calculation of the Spin Stiffness in the J1—J2 Heisenberg Antiferromagnet 
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We calculate the spin stiffness p s for the frustrated spin-| Heisenberg antiferromagnet on a 
square lattice by exact diagonalizations on finite clusters of up to 36 sites followed by extrapolations 
to the thermodynamic limit. For the non-frustrated case, we find that p s = (0.183 ± 0.003)Ji, in 
excellent agreement with the best results obtained by other means. Turning on frustration, the 
extrapolated stiffness vanishes for 0.4 < J2/J1 ^ 0-6. In this intermediate region, the finite-size 
scaling works poorly - an additional sign that their is neither Neel nor collinear magnetic order. 
Using a hydrodynamic relation, and previous results for the transverse susceptibility, we also estimate 
the spin-wave velocity in the Neel-ordered region. 
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The question of the existence of long-range magnetic 
order (LRMO) in systems with frustrated interactions 
and strong (quantum or thermal) fluctuations is often 
difficult to decide. The traditional way of answering 
this question is by calculating magnetic order parame- 
ters. An alternative way is to consider the spin stiffness 
p s , which is non-zero in a LRMO state. The stiffness has 
the advantage of being unbiased with respect to the order 
parameter, and constitutes, together with the spin-wave 
velocity, the fundamental parameter that determines the 
low-energy dynamics of magnetic systems [1]. It is there- 
fore of importance to find accurate values for p s . 

The spin stiffness measures the energy cost to intro- 
duce a twist 6 of the direction of spin between every pair 
of neighboring rows, 
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where Eo(9) is the ground-state energy as a function of 
the imposed twist, and N is the number of sites. In 
the thermodynamic limit, a positive value of p s means 
that LRMO persists in the system, while a zero value 
reveals that there is no LRMO, as is the case in a spin- 
liquid. When looking at a finite system, things are more 
complicated. Here the stiffness is only zero at distinct 
points, and is positive or negative on the intervals in be- 
tween. A negative value says that the system is unsta- 
ble to a change in the boundary conditions, suggesting 
that the true ground state of the model in the thermo- 
dynamic limit is incommensurate with the structure of 
the finite cluster being used. A positive value reveals 
a stable ground state, and can sometimes be used with 
finite-size scaling to extract the behavior of the stiffness 
in the thermodynamic limit. This is in particular the 
case in the Neel and collinear regions. 

The spin stiffness for the unfrustrated spin--| Heisen- 
berg model on a 2D square lattice has been calculated 
directly by series expansion [2], p s = (0.18 ± 0.01)Ji, 
and by second-order spin- wave theory (SSWT) [3], p s = 
(0.181 ± 0.001) J\. Furthermore, the spin-wave velocity c 
and the transverse susceptibility x± have also been calcu- 



lated in SSWT, and since the ensemble of values fulfill the 
hydrodynamic relation [1] p s = c 2 \± to a good approxi- 
mation, there is strong evidence for the accuracy of the 
SSWT values [3]. However, a previous attempt to extract 
the value of p s from exact diagonalizations (ED) yielded 
p s = 0.125Ji [4], which is far away from the other results. 
This is not too bothering regarding that the ED value of 
p s (and c) was obtained from the correction terms in the 
finite-size scaling analysis and as such looses accuracy 
due to cancelations, and is further influenced by higher- 
order corrections which are not known. To obtain more 
accurate values of the spin stiffness, we here set out to 
calculate the spin stiffness directly by using EDs to evalu- 
ate p s as a correlation function. In contrast to two recent 
works which have employed ED and finite twists on the 
square and triangular lattices [5,6], our method preserves 
more symmetries, and we can treat clusters of up to 36 
sites. 

By performing a careful finite-size extrapolation we ar- 
rive at a value of the stiffness in the non-frustrated case, 
p s = (0.183 ± 0.003)Ji, in excellent agreement with the 
SSWT and series-expansion results. In the case of frus- 
trating interactions, things are more complicated. In a 
previous ED study [4], the order-parameter was found to 
vanish in the region 0.34 < J2/J1 ~ 0.68, and one of our 
aims was to find out whether a direct calculation of the 
stiffness would corroborate this result. Our results sug- 
gest that the stiffness vanishes for 0.4 < J2/J1 ^ 0.6, but 
there is also a tendency of the stiffness to blow up in the 
region J 2 / Ji ^ \ ■ A similar tendency is found in a first- 
order SWT (FSWT). In the latter case, this burst is a 
signature of the breakdown of SWT as J 2 / J\ approaches 
the classical transition point J2/J1 = \- 

We start with the general Heisenberg Hamiltonian 
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where the sum goes over all pairs of sites and in- 

troduce a local rotation at site i by 0i around the z- 
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that Sj ot is unchanged. A Mac-Laurin expansion around 
9ij = 9i — 9j =0 gives to order [5] 
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^Jij(S^"S- — S i S~j~) is the z-component of 



the spin-current operator, and T{j = ^J{j(SfSj + S~S+) 
is the "spin-kinetic-energy" operator. To obtain the spin 
stiffness, a uniform twist 6 is introduced between each 
pair of adjacent rows, i.e., Oij = #[(r; — rj) ■ y], and to 
second order in 6 one has (H(6)) = (Ho) + ^N6 2 p s . This 
gives a direct expression for p s , which for the J\-J'i model 
[J 8J - = ,J 1 (J 2 ) for nearest (next-nearest) neighbors] reads 
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where the JY terms comes from second order perturba- 
tion theory, and where P = I — |0)(0| is the projec- 
tion operator projecting on the space orthogonal to the 
ground state. Note that the J 2 term has two terms per 
site and that the expectation values are evaluated in the 
wow-twisted space. The stiffness is now expressed as a 
sum of a "kinetic-energy" term TY, which is easy to cal- 
culate, and a spin-current spin-current correlation func- 
tion JY, which needs some computational efforts to be 
evaluated. 

To calculate JY, we use a continued-fraction expan- 
sion [7] where we repeatedly apply the Hamiltonian on 

the spin-current state |/o) = Pojy S ^\0), which is antisym- 
metric under spin inversion and under reflection on the 
a;-axis. The loose of diagonal reflection symmetry im- 
plies a doubling of the size of the Hilbert space, which 
for the 36-site cluster is now ~ 3 x 10 7 . The expansion 
normally converges very quickly and p s is obtained with 
five significant digits after five to ten iterations. 

As a test of our method, we first considered the fer- 
romagnetic model, where both J\ and J 2 are negative. 
The ferromagnetic state with S* ot = is the symmet- 
ric superposition of all S* ot = states. The transverse 
correlations are easily obtained as (i(S^~S~ + S~S^~)) = 
^N/(N — 1) and, for periodic boundary conditions, the 
JY term is identically zero. The order parameter lies 
in the z = plane, and one is really measuring the 
(transverse) spin stiffness (compare the AF case below), 
p s = ^(Ji + 2J2)N/(N — 1). This result is exactly repro- 
duced in our exact diagonalizations. 

The antiferromagnetic case differs from the FM case 
both by the necessity to consider the spin-current term 
and by the ground state being rotationally invariant. The 
latter fact means that the twists are not orthogonal to 
the order parameter, but instead we calculate the rota- 
tional average of the stiffness. Since the stiffness for a 
twist around the Neel (or the collinear) order parameter 
is zero, we have to multiply our result by a factor | to 
arrive at the ordinary transverse stiffness. Let us first 



consider the unfrustrated case. 

To extract the values of thermodynamic quantities 
from finite-size calculations it is of crucial importance 
to have good knowledge about the scaling behavior of 
the quantities of interest. A great deal of information 
can be obtained from studying how the spin-wave theory 
behaves under scaling, or from the finite-size analysis of 
the non-linear a model [8]. The FSWT expression for 
the stiffness [3] can be written as 
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where Eo is the LSWT ground-state energy and is re- 
lated to the LSWT dispersion relation by lo^ = ASJe^. 
By looking at the fc-sums involved, one finds that the cor- 
rection to the ground state energy per site Eo/N scales as 
N~ 3 I 2 and that the correction to the second term scales 
as TV -1 / 2 . Using the rotational invariance of the ground 
state, we can further rewrite the ED expression (4) as 



The physical content of the first term is thus exactly the 
same in the both cases, and it is known that the cor- 
rection to Eo/N goes as TV -3 / 2 also in the ED case [8]. 
It is therefore wise to use the same scaling as in SWT 
also for the JY term, 3Y N - JYoo oc N' 1 / 2 . With these 
scaling laws we can extrapolate the TY and JY terms 
separately, and then finally obtain the stiffness in the 
thermodynamic limit as 

Ps,oo — TYqq -\- JYqq . (9) 

As was noted in Ref. [4], the extrapolated value is sen- 
sitive to which set of cluster sizes one uses. In Tab. I, 
the results for the different clusters are presented and in 
Tab. II, the results of the various extrapolations are pre- 
sented together with error bounds coming from a x 2 -fit 
of the values to a straight line. As seen in Tab. II, the set 
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of clusters with {18, 20, 32, 36} sites gives the best result 
in the non-frustrated case. When turning on J 2 we are 
in a much less understood regime. Semi-classically, there 
is a sharp transition at J2/J1 = \, from a Neel state to 
a collinear state. However, going to S = there may 
well be a widening of the transition region and a region 
with a spin-liquid ground state may open up. Indeed, 
the earlier finite-size studies suggested that the Neel and 
collinear states are separated by an intermediate region 
for 0.34 < J 2 /Ji < 0.68 [4]. On the other hand, be- 
sides a number of works which have also found a reduced 
Neel stability, the large-5 1 studies using Schwinger-boson 
mean-field theory [9] or a self-consistent spin-wave theory 
[10] show an increased Neel stability with respect to the 
classical case. Since these methods are only trust-worthy 
for large values of S, the discrepancy for S = ^ is not 
necessarily significant. It is also not surprising that a self- 
consistent mean-field calculation of p s yielded a stiffness 
which does not vanish until J2/J1 ~ 0.6 [11]. 




0.0 



-0.2 



0.0 



0.2 



0.4 

J 2 /J, 



0.6 



0.8 



1.0 



FIG. 1. The stiffness p s for the various clusters being used. 
The 18-site cluster shows a negative stiffness for big J2, and 
the 20-site cluster has a change in the ground-state symmetry 
around J2 / Ji = 0.58. 

A good test of our numerical program is to consider 
the limit J2/J1 = 00, or J\ = 0, J2 = 1. Here, the two 
sublattices decouple and the energy (stiffness) should be 
twice the energy (stiffness) of the subclusters. This is 
indeed exactly what we obtain. As J2/J1 increases, we 
thus expect to see a decrease in the stiffness followed by 
an increase as the two sublattices become individually 
ordered. The minimum should be somewhere not too far 
from the classical break point J2/J1 = \- For the 18-site 
cluster the stiffness should go negative for large J 2 / Ji 
because that cluster is not compatible with the structure 
of two antiferromagnetic sublattices. These observations 
agree with the results for the finite clusters presented in 
Fig. 1. 

Unfortunately, the individual properties of the clusters 
result in rather strong peculiarities. In the 20 and 36-site 
cases, there is a change in the symmetry of the ground 
state as the two sublattices become individually antifer- 
romagnetically ordered. In the 20-site case, this causes 



an abrupt jump in the stiffness, while in the 36-site case 
the transition is very smooth. 
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FIG. 2. The finite-size data for TY and JY for J 2 /Ji = 
(crosses), 0.2 (diamonds), and 0.4 (circles) together with fits 
for J2 / Ji = 0.0 and 0.2. For small frustration, J2 / Ji < 0.3, 
the scaling law is well satisfied, but in the intermediate region 
the results do not line up. 

Given the strong individual differences in Fig. 1, it is 
not evident how to extrapolate to N = 00 for the various 
degrees of frustration. In Fig. 2, we show the actual data 
which we try to fit with our scaling laws, for J2/J1 = 0, 
0.2, and 0.4. In the region 0.3 < J2/J1 ~ 0.6, the results 
for JY do not line up and the extrapolation to N = 00 
is unreliable. In Fig. 3, we show the results of extrapo- 
lations using a few different sets of clusters. In the in- 
termediate region our results are scattered. The FSWT 
result is obtained by generalizing Eq. (7). 




FIG. 3. The extrapolated value of the stiffness for 
some choices of clusters together with earlier ED results 
(ED+NL(tM) and FSWT. 

By excluding the 20-site cluster, the results suggest 
a vanishing LRMO in the region 0.4 < J2/J1 ~ 0.6 in 
rather good agreement with the previous ED results [4] 
(where the stiffness vanished at the same point as the or- 
der parameter). The extrapolation from the {20, 32, 36}- 
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site clusters follows the FSWT result closely, but the 
derivative of the latter diverges as J 2 / J\ — > \ and we 
consider this similarity to be fortuitous. If one really 
were in the Neel regime all the way to J1/J2 ^ \, the cou- 
pling constant in the non-linear a model, g oc c/p s , would 
be roughly constant over the entire region and there is 
no reason why the finite-size scaling should cease to be 
valid. This is however the case as seen in Fig. 2, and we 
conclude that the intermediate region has neither Neel 
nor collinear order, and that a first-order transition from 
Neel to collinear order as suggested in Refs. [9] and [10] 
is inconsistent with this result. 

Since we consider our result from the {16,32,36}- 
cluster extrapolation to be good, we can combine it with 
the previous ED results for the transverse susceptibility 
X±_ [4] to obtain the spin-wave velocity c from the hy- 
drodynamic relation c = \J p s /x± • The result is shown 
in Fig. 4. The result i s in fair agreement with LSWT, 
c = Ji \/2(l — 2 J2I Ji), but close to the phase bound- 
ary the result may not be trusted since the susceptibility 
and the stiffness do not vanish at the same point. In the 
non-frustrated case, our best value, p s = 0.183Ji, yields 
c = 1.67Ji in excellent agreement with the SSWT result 
[12,3] c= 1.664Ji. 

Bonca et al. [5] have reported results for p s for the 16 
and 20-site clusters. Their results differ from ours due 
to a number of lapses on their side. First of all, they 
did not include the J2 terms in Eqs. (5,6). Secondly, 
they missed the factor |, which compensates for the ro- 
tational symmetry of the ground state, and finally they 
did not use the proper power laws in their extrapolation 
to the thermodynamic limit. 

It would be of great interest to extract some precise 
signature of the ground state in the intermediate region. 
This is however not possible from the spin stiffness. Even 
a spin liquid may have a finite stiffness for a finite sys- 
tem and in the region where the finite-size scaling does 
not work, we can only exclude Neel and collinear LRO. 
Our results strongly suggest the existence of an uncon- 
ventional ground state in a wide intermediate region, but 
its nature has to be revealed by a more detailed exami- 
nation of the correlation functions. 




FIG. 4. The spin-wave velocity obtained by using the hy- 
drodynamic relation. As a comparison, the LSWT results are 
shown. 
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TABLE I. Values of p s , TY, and 1Y for finite clusters and 
no frustration. 



TY 



1Y 



Ps 



16 
18 
20 
32 
36 



0.35089 
0.34699 
0.34540 
0.34009 
0.33943 



-0.07248 
-0.08054 
-0.08433 
-0.09938 
-0.10084 



0.27841 
0.26646 
0.26107 
0.24071 
0.23859 



TABLE II. Extrapolated values for TY, IY, and p s for 
J2 = 0, with the uncertainty in the last digit given in paren- 
theses. 



Cluster sets 



TY 



JY 



16,18,20,32,36 

16,20,32,36 

16,32,36 

18,20,32,36 

20,32,36 



0.3345(7) 
0.3344(6) 
0.3344(2) 
0.3352(1) 
0.3351(2) 



-0.157(5) 
-0.159(6) 
-0.160(5) 
-0.152(3) 
-0.152(5) 



0.177(6) 
0.176(7) 
0.174(5) 
0.183(3) 
0.184(5) 
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